# generate graphs of coefficients #

# estimate effects (raw) #

# run regressions #

coef = rep(0,length(depvars)) # vector to store coefficients #

se = rep(0,length(depvars)) # vector to store standard errors #

# run regressions #

for (i in 1:length(depvars)) {
  
  f = as.formula(paste(depvars[i],treatment, sep = " ~ "))
  model = lm_robust(f, data = agro, se_type = "HC1")
  coef[i] = model$coefficients[2]
  se[i] = sqrt(diag(vcov(model))[2])
  
}

# store coefficients and standard errors as text #

text = rep("",length(depvars)) # vector to store text #
beta='\u03b2'

for (i in 1:length(text)) {
  
  text[i] = paste0(beta," (se) = ", 
                   formatC(coef[i], format = 'f', flag='0', digits = 2)," (",
                   formatC(se[i], format = 'f', flag='0', digits = 2),")")
  
}

# estimate effects (normalized) #

agro = agro %>% mutate_at(depvars,scale)

# data to store results #

spec =  as.data.frame(rep(2:1, each = length(depvars)))
number = as.data.frame(c(1:length(depvars)))
colnames(spec) = "model"
colnames(number) = "depvar"

# specification 1 #

coef = rep(0,length(depvars)) # vector to store coefficients #
se = rep(0,length(depvars)) # vector to store standard errors #
for (i in 1:length(depvars)) {
  f = as.formula(paste(depvars[i],treatment, sep = " ~ "))
  model = lm_robust(f, data = agro, se_type = "HC1")
  coef[i] = model$coefficients[2]
  se[i] = sqrt(diag(vcov(model))[2])
}
results1 = cbind(number,coef,se)
results1 = results1 %>% mutate(icmax = coef+1.96*se, icmin = coef-1.96*se)

# specification 2 #

coef = rep(0,length(depvars)) # vector to store coefficients #
se = rep(0,length(depvars)) # vector to store standard errors #
for (i in 1:length(depvars)) {
  f = as.formula(paste0(depvars[i]," ~ ",treatment,c(" + yob + sex + res + main + coop + area")))
  model = lm_robust(f, data = agro, se_type = "HC1")
  coef[i] = model$coefficients[2]
  se[i] = sqrt(diag(vcov(model))[2])
}
results2 = cbind(number,coef,se)
results2 = results2 %>% mutate(icmax = coef+1.96*se, icmin = coef-1.96*se)

# data #

results = rbind(results1,results2)
results = cbind(spec,results)
results$text = text

# plot coefficients #

graph = ggplot(results, aes(x = depvar, y = coef, colour = as.factor(model))) + 
  geom_hline(yintercept = 0, colour = "#323232", lty = 2) + 
  geom_point(size = 2, position = position_dodge(width = 1/2)) +
  geom_linerange(aes(ymin = icmin, ymax = icmax), size = 1, position = position_dodge(width = 1/2)) +
  theme_linedraw() + 
  xlab("") + ylab("Standardized Effects") + 
  theme(plot.title = element_text(size=14, family="sans", hjust = 0.5, face = "bold"),
        axis.title.x = element_text(size=13, family="sans", face = "bold", colour = "#323232"), 
        axis.title.y = element_text(size=13, family="sans", face = "bold", colour = "#323232"),
        axis.text.x = element_text(size=11, family="sans", colour = "#323232"), 
        axis.text.y = element_text(size=12, family="sans", colour = "#323232"),
        #panel.grid.major =  element_blank(), 
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "#323232", size = 0.6), 
        axis.ticks = element_line(colour = "#323232", size = 1),
        plot.margin = margin(1, 0.5, 0.5, 0.5, "cm"),
        legend.position = "bottom", legend.title = element_blank(), 
        legend.text = element_text(size=12)
  ) +
  scale_x_continuous(breaks=1:8, labels = c("log(Expenditures)",
                                            "Pesticide Use (0/1)", "Tractor Use (0/1)", 
                                            "Conserv. Practices (#)", "Manag. Practices (#)", 
                                            "Rotational Grazing (0/1)", "Pasture Restoration (%)", 
                                            "Technical Assistance (0/1)")) + 
  scale_y_continuous(breaks = c(-0.4,0,0.4,0.8,1.2,1.6,2.0), limits = c(-0.4,2.4)) + coord_flip() +
  geom_text(aes(x=depvar, y=2, 
                label = text),
            size = 4, hjust=0.35, vjust=-0.14, colour = "#323232", family="sans") + 
  scale_color_manual(values=c("#d7191c","#fdae61"), 
                     breaks = c("2","1"),
                     labels = c("Bivariate","Demographic Controls"),
                     name="")

